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We present a detailed study of the Equilibriumlike invaded cluster algorithm (EIC), recently 
proposed as an extension of the invaded cluster (IC) algorithm, designed to drive the system to 
criticality while still preserving the equilibrium ensemble. We perform extensive simulations on two 
special cases of the Potts model and examine the precision of critical exponents by including the 
leading corrections. We show that both thermal and magnetic critical exponents can be obtained 
with high accuracy compared to the best available results. The choice of the auxiliary parameters 
of the algorithm is discussed in context of dynamical properties. We also discuss the relation to the 
Li-Sokal bound for the dynamical exponent z. 
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I. INTRODUCTION 

As an important tool in the study of phase transitions 
Monte Carlo simulations passed through significant ad- 
vances in the last decades and number of new algorithms 
and improvements were proposed in various directions, 
such as cluster algorithms by Swendsen and Wang [l[ 
or Wolff [H multicanonical algorithm 0], Wang-Landau 
algorithm jj]. An interesting direction undertaken by 
Machta et al. [|[ was to construct an algorithm which 
would, in a self organized way, drive the system to criti- 
cal temperature without prior knowledge of it. It was ap- 
plied to various classical models from discrete to continu- 
ous ones including some more complicated effects such as 
frustration @, quasi periodic ordering 0, or tricritical 
points . Their method, based on the invasion percola- 
tion Q and the cluster algorithm, appeared to converge 
much faster than in the Swendsen- Wang (SW) algorithm, 
with a dynamical exponent z close to zero. In the same 
time it opened a number of issues [l(J,[ll| related to the 
fact that the underlying nonequilibrium procedure does 
not generate an equilibrium ensemble and consequently 
is not able to reproduce correct critical exponents when a 
finite size scaling is applied. There were also some other 
attempts to design a cluster algorithm that would self- 
regulate to criticality [IH, EU . 

Recently we proposed a modification of the IC algo- 
rithm [l4j aimed to reestablish the equilibrium in the self- 
regulating procedure of the IC algorithm. The principal 
idea is to impose a simple constraint on the temperature 
uncertainty characteristic to the IC algorithm, reducing 
it to the limits compatible with the equilibrium distribu- 
tion. As shown in our preliminary work (l4j . when this 
single constraint is applied to the IC algorithm, the cor- 
rect scaling properties of thermodynamic obscrvablcs are 
recovered, while the algorithm still retains the property 



of self driving to criticality. 

The aim of the present paper is twofold. First, it of- 
fers a more intensive and comprehensive numerical cal- 
culations, including the leading convergence exponent, in 
order to examine the precision of the results for critical 
exponents and temperature that the EIC algorithm can 
provide. Second is to extend the study to the dynami- 
cal properties of the algorithm, in particular to the au- 
tocorrelation functions, which indicate deviations from 
the equilibrium distribution providing the criteria for 
the choice of auxiliary parameters of the method, which 
would not produce systematic errors. 

The paper is organized as follows. In Sect. |TT]we in- 
troduce the model and notations and describe in detail 
the principles of EIC algorithm. In Sect. Illll is presented 
the analysis of critical exponents and temperature in- 
cluding the leading correction to scaling. We also add 
a brief comment on the application of this algorithm to 
the first-order phase transition. In Sect. [IV] we study 
the dynamical properties of EIC algorithm through the 
autocorrelations of energy and order parameter with a 
special emphasis on the impact of choice of the two aux- 
iliary parameters. Sect. IVl contains the conclusion. 



II. ALGORITHM 

We consider the Potts model [IH [l6( defined by the 
Hamiltonian 

H = - J E (w,-i)> (!) 

<i,j> 

where Oj denotes the g-state Potts variable at the lattice 
site i and the summation runs over the nearest neigh- 
bors. As shown by Fortuin and Kasteleyn (FK) [T3|, the 
partition function of the model ([I]) can be written as 
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where 



p=l 



(3) 



and /3 is the Boltzmann factor. By expanding the bino- 
mial product into sum over graphs and integrating over 
the spin degrees of freedom it reduces to the random- 
cluster (RC) model 
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which can be understood as a generalized bond percola- 
tion, where p is the bond probability. The summation 
in Eq. Q is taken over the set T of all the graphs on 
the lattice, each graph 7 representing one possible bond 
configuration. B is the number of the lattice edges, 6(7) 
denotes the number of bonds, and 0(7) is the number of 
connected components (FK clusters) in the graph 7. 

FK expansion of the Potts model partition function is 
also used as a basis for construction of cluster algorithms 
for numerical simulations [l8|. To this purpose the 
r.h.s. of the Eq. ^ can be written in alternative way by 
introducing new degrees of freedom m t j — 0, 1, represent- 
ing the absence or presence of bonds in the expansion of 
the binary product in Eq. ([2|). In such a way the Boltz- 
mann weight in ([2]) is decomposed into two joint distri- 
butions: over spin ({cr}) or bond ({«}) configurations 
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As it may be seen from Eqs @ and for a given 
bond configuration, only the spin configurations where 
the spins belonging to the same cluster are in the same 
state give a nonzero contribution to the partition func- 
tion, while entire clusters may be flipped without an en- 
ergy cost. On the other hand, for a given spin configura- 
tion ({cr}) the Boltzmann weight reduces to a binomial 
distribution of bond occupancy 



({a}) = L(l- P f-^.J2 h S ) P b -(l- P r- b , (6) 



where n ss denotes the total number of "satisfied" edges 
(surrounded by spins in the same state). The standard 
cluster algorithms such as SW [7J consist in alternate up- 
dates of bonds and spins. In the first part of the MC 
step, taken a given spin configuration {cr}, the FK clus- 
ters are formed by putting bonds, with probability p, 
between neighbors in the same state. In the second part 
each cluster is flipped at random producing a new con- 
figuration of spins. The bonds are then erased and one 
proceeds to the next MC step. In such a way the equi- 
librium ensemble is generated. 

The IC algorithm Q consists of the same two steps, 
but the bond update is altered in order to accomplish 
the self driving to criticality. Given the configuration of 
spins, the bonds arc added with the probability 1, but 



only until the percolation is achieved. This is interpreted 
as a signature of a finite-size critical point. As soon as the 
percolation is reached the FK clusters arc randomized. 
The fraction of "satisfied" edges that has been populated 

KW}) 



P{a} 
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n ss {{(r}) ' 

is recorded after each MC step. b({a}) denotes the num- 
ber of added bonds and n ss ({a}) the number of neigh- 
boring pairs in a same state for a given spin configuration 
{cr}. The average over simulations P{ a } is identified with 
the quantity p defined in (|3|), yielding a critical temper- 
ature estimate. It is easy to understand why such a pro- 
cedure is self-driving to criticality. Namely if the given 
configuration of spins is typical for a temperature much 
higher than T c (with random and poorly correlated spins) 
almost all the edges with satisfied neighbors will need to 
be populated before the percolation has been achieved. 
The n ss in the next iteration can only be larger - the 
system will be driven towards configurations with lower 
energies, corresponding to lower temperatures. On the 
other hand, when the configuration is typical of lower 
temperatures, with n ss considerably larger than the one 
in criticality, small amount of added bonds will suffice to 
achieve percolation, so that, after randomization of clus- 
ters n ss will be diminished driving the system to higher 
temperature configurations. In the same time one can 
notice here a back-draw of this self driving procedure. In 
both cases, the described configurations that have been 
retained for statistics result from the tails of a binomial 
distribution ([6]). Problem is that such configurations are 
typical during entire IC simulation. Namely, in spite 
of the fact that the system is driven to criticality, the 
described oscillations remain significant throughout the 
simulations, so that it generates distribution of P{ a } that 
is much wider than it would result from the binomial dis- 
tribution of ([6]). Consequently, the resulting ensemble is 
not a canonical one and not in equilibrium, as observed 
already in jToj . 

Our EIC algorithm approach follows the same IC pro- 
cedure but restrains the excessive fluctuations in P{ a } 
in order to regain the equilibrium fluctuations. The 
main idea is to constrain the width of distribution of 
P{ a } obtained from the simulations by the relation ([7} 
to the width which is compatible with the distribution 
of b{ a y in (|6|) and which has the same scaling in L. 
The width of the binomial distribution in Eq. ([6]) gives 
var{jp{ a }) oc y/p(l — p)n ss d ^' 2 . Since n ss corresponds to 
the energy of a given configuration and is oc L d , we im- 
pose the constraint on P{ a }, which allows it to vary only 
within limited range of the width v oc L~ d / 2 . 

To include this constraint into algorithm, we modified 
the stopping rule by introducing two auxiliary parame- 
ters to the IC algorithm, v and N a . Let us review this 
procedure briefly and discuss the role of these free pa- 
rameters. 

The MC iterations are grouped in intervals of N a MC 
steps. In every interval i of N a steps the bonds are placed 
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on the lattice between satisfied neighbors with the inten- 
tion to construct a percolating cluster. However, the sys- 
tem is not left to percolate with just any value of P{ a } 
as in the IC algorithm. Any accepted configuration is 
required to lie inside the interval: — v < ps a \ < 

p~i_i + v, where p i denotes the average of P{ a } over the 
steps. Here v is the free parameter set 



i-th block of N„. 
to be: 



v = v- L- d ' 2 



(8) 



and v is a constant whose choice will be discussed later. 
If the system percolates before the lower bound of P{ a } is 
reached, the bonds are still added until the lower bound is 
attained. If the upper bound of P{ a ] is reached before the 
percolation is attained, the process is stopped without 
requiring the system to percolate. After every block of 
N a MC steps, the new average, p 4 , is calculated to be 
used as a reference value in the next N a steps. Obviously, 
in this way the information on the value oi p i _ 1 in the 
preceding N a MC steps is propagated to the consecutive 
interval i and will have an effect on the autocorrelation 
functions (to be discussed in the Sect. IIV|) . 

The starting value p can be obtained by the uncon- 
strained IC algorithm. One can also define a faster equi- 
libration of the IC algorithm by gradually reducing the 
value of parameter v in (jHJ during some initial set of in- 
tervals of N a steps. These first few blocks of N a steps 
are discarded from statistics and the configurations are 
recorded after the steady state was reached. 

Distribution of P{ a } resulting from the above proce- 
dure has the required width proportional to v = v ■ L~i 
under the assumption that N a is large enough so that the 
fluctuations of P{ a } within the interval are dominant over 
the fluctuations of p i between the intervals. In this case 
the fluctuations of the mean value p i reduce to the actual 
temperature fluctuations, producing the width of order 
L~ x l v < L~ d / 2 if the heat capacity exponent a > 0, or 
of order L~ d l 2 if a < 0. If N a would be small enough 
for the fluctuations of p.^ between the intervals to take 
over, a crossover to the IC regime would be seen. Later 
we will demonstrate that there exists a clear indication 
when this occurs. 

In Fig. [T] we give an illustration that the canonical 
form of energy distribution is recovered within the EIC 
algoritm. As noted by Machta etal. 0, L d ■ var(e) for 
2D Ising obtained by the IC algorithm scales approxi- 
mately linearly with the system size L, while it should 
scale as InL if the ensemble was canonical. To examine 
the same point within the EIC algorithm, we present in 
Fig. [I] an overlap fit of the energy distributions resulting 
from EIC algorithm and compare it with the distribu- 
tions obtained by IC algorithm with equal statistics, for 
several lattice sizes in the case of q = 3 Potts in 2D. The 
distributions are rescaled with the exact value of the ex- 
ponent describing the width of the energy distribution in 
the equilibrium ensemble, X = | — = 0.8. It is im- 
mediately seen that the results for different L, obtained 
by EIC, collapse very well, while the ones obtained by 
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Figure 1: (Color online) Overlap fits of energy distributions 
w(e) for q — 3 Potts in 2D rescaled with the exact exponent 
X = 0.8. In the inset are shown ln-ln plots of L d ■ var(e) 
vs. L. The statistic is 10 6 MCS for both EIC and IC cases, 
and the parameters of the EIC algorithm were chosen to be 
N a = 100 and v = 



IC do not. To quantify that fact, we show in the inset of 
Fig. Q]the ln-ln plots of L d ■ var(e) vs. L, the slopes of 
which in the canonical ensemble represent the exponent 
otjv (exactly equal to 0.4 in this case). For the lattice 
sizes considered we obtain rj 0.38 and ps 1.26 for EIC 
and IC cases respectively (in terms of the width scaling 
Xeic ~ 0.81 and Xjc ~ 0.37). This indicates that the 
ensemble sampled by the EIC algorithm is indeed canon- 
ical, contrary to the one produced by the standard IC 
algorithm, which appears to be much wider and scales 
with different exponent, characterizing the algorithm it- 
self [To|. As it shall be seen in the next section, the EIC 
algorithm can be used to obtain very accurate values for 
critical exponents, not only the magnetic one related to 
criticality, but also the thermal critical exponent related 
to the approach to criticality and implying the validity 
of fluctuation-dissipation theorem. 



III. CRITICAL BEHAVIOR 

Derivation of critical exponents and temperature by 
EIC algorithm was already illustrated by applying it to 
several cases of the Potts model in our preliminary study. 
In present work, we selected only two cases both ex- 
hibiting 2nd order phase transition, one in two and one 
in three dimensions: case q = 3 in 2D and Ising case 
(q = 2) in 3D. We conducted extensive numerical simu- 
lations increased by an order of magnitude with respect 
to the preliminary study 



14j . This allowed us to calcu- 



late the leading power-law corrections as well and analyze 
the precision of the results with particular attention to 
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possible systematic errors due to the choice of the auxil- 
iary parameters. Unlike our previous study, we avoided 
to use any data exactly known in advance and applied 
consistently multi parameter fits, which give critical tem- 
perature and critical exponents simultaneously whenever 
needed. 

The simulations were performed on lattices with pe- 
riodic boundary conditions and sizes up to L = 96 and 
L = 1024 in the 3D and 2D case respectively. The statis- 
tics varied from 15 • 10 6 iterations for smaller to 7 • 10 6 
iterations for the largest lattice sizes in the 3D case and 
from 15 • 10 6 to 5 • 10 6 in the 2D case. The auxiliary 
parameters were set to v = 1/10 and N a = 100 in both 
cases. The percolation was established by the topolog- 
ical rule, i.e. by the condition that the infinite cluster 
wraps around the lattice. The time necessary to reach 
the equilibrium state never exceeded the first 3000 MC 
steps with the above choice of auxiliary parameters N a 
and v. 

The running time per MC step was approximately the 
same as for the original IC algorithm. For illustration, a 
run of 10 5 MC steps for the 2D Ising model on L = 64 
lattice requires approximately 500s on the AMD Opteron 
240 processor (1.4GHz). 



A. The magnetic critical exponent yn 

The magnetic critical exponent was evaluated indepen- 
dently from three different quantities calculated at criti- 
cality, the average mass of the largest cluster, the order 
parameter and the susceptibility. 

The average mass of the largest cluster s max is defined 
as the number of spins in the largest cluster averaged 
over all the configurations. Its scaling yields the critical 
exponent directly 
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A four-parameter fit of our finite-size results (see Fig. 
[2]) to the form ([9|) produced rather accurate estimate of 
the exponent yh and showed that the leading correction 
is of a power law form. Results for the exponents are: 
y h = 2.4815(5) with u t = -2.2(1) and y h = 1.8661(7) 
with u>i = —1.1(1) for 3D Ising and 2D q=3 Potts model 
respectively. 

The order parameter m scales at T c according to 



\P=Pc 



(10) 



where (3/v is related to yh by ft/v = d — y%. It was 
calculated by using the standard definition of the order 
parameter for the Potts model 



m = mi, 
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which consists in taking the most populated among the q 
Potts states in each configuration. By a four-parameter 
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Figure 2: (Color online) Largest cluster mass ~s m ax versus the 
system size. 



fit to the form (JTOJ) we found the exponents fi/v = 
0.5180(2) with uj 2 = -1.9(1) and fi/v = 0.1345(5) with 
uj 2 = —2.2(1) for 3D Ising and q=3 case in 2D respec- 
tively. 

The fluctuations of the magnetization in equilibrium 
are related to the susceptibility by the relation 
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m 2 — m 2 
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which assumes the validity of the fluctuation-dissipation 
theorem, not fulfilled for the standard IC algorithm. For 
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Figure 3: (Color online) Order parameter fluctuations plotted 
versus the system size L. 

a finite system, at the quasicritical point T c l, the suscep- 
tibility exhibits maximum, which scales with a power-law 
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and is related to the magnetic exponent by 7/1/ = 2 • 
Uh — d. The log-log fit of the results (sec Fig [3|) gives 
7/V = 1.98(1) and -y/v = 1.74(1) for 3D Ising and q=3 
case in 2D respectively. We find the values of j/v less 
accurate than the estimates for y^ or ft/v. The reason 
may be attributed to the fact that, for a finite system, 
the critical point identified by the percolation through a 
specific choice of the stopping rule, should not necessarily 
coincide with the quasi-critical temperature defined by 
the maximum of the susceptibility. 

Another quantity of interest, related to the magneti- 
zation is the Binder fourth order cumulant fl9l ]. defined 
as 



results arc shown for 3D Ising and 2D q = 3 cases. 
The values can be fitted quite accurately to a three- 
parameter power law form: U±(L) = [/4(a) + bjj.i ■ L" 3 . 
The constants 1/4(0,) are found to be /74(a) = 0.412(4) 
and U A (a) = 0.523(2) for 3D Ising and q=3 Potts in 2D 
respectively. The respective convergence exponent CJ3 is 
determined to be 013 = —0.75(3) and W3 = —0.91(4). A 
value Ui(L 00) « 0.40 which is similar to the value 
obtained in the present work for 3D Ising case can be 
found in Ref. [21| which was obtained by the extrapo- 
lation of th finite-size values U^L) taken at the point 
where dU Qj L ^ attains its maximum. 



U 4 = 1 
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One of its main benefits within numerical calculations 
is to provide a good estimation of the critical temper- 
ature, which in the present case, where the algorithm 
drives itself to the critical point is not the main issue. At 
criticality it may be related to the one of the universal 
amplitude ratios, and we found interesting to examine 
this quantity and its size convergence within the present 
algorithm. 

It is important to notice, that in the present approach, 
U4 is calculated at the finite size critical point defined by 
the topological rule, hence t v L = a, where a is some non 
universal constant factor. The L dependence of Binder 
cumulant at this finite size critical point is 



U 4 (L, t, u) = Ui(T v L) + b lU L y ' + 



(15) 



Consequently, the constant term of the Binder cumulant 
in the thermodynamic limit depends on the way of ap- 
proaching to the critical point. Other authors [20( con- 
sidered instead the universal value 1/4(0). In Fig [4] the 
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Figure 4: (Color online) EIC data for Binder cumulant Un(L) 
versus 1/L. 



B. The thermal critical exponent y t and critical 
temperature 

Two quantities related to the thermal critical exponent 
y T — 1/v were examined. The first one is p c (L), related 
to the quasicritical temperature, with the leading size 
dependence 



p c {L) = p c (L oo) 



(16) 



It is obtained from the average over P{ a } defined in the 
Eq. ([7]). Data for p c (L) with a three-parameter fit are 
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Figure 5: (Color online) Data for p c (L) versus 1/L. Dotted 
lines describe three parameter fits and crosses denote exact, 
or best known values. 

illustrated in Figure [5l The fit appears to be more suit- 
able for finding the p c {L oo), giving the critical point 
with the precision up to the sixth digit, yielding p c (L — > 
oo) = 0.358097(1) and p c (L oo) = 0.633975(1) for 
3D Ising and q=3 case in 2D respectively. An estimate 
of the exponent y t by the same procedure appears, how- 
ever, to be by an order of magnitude less precise than by 
the logarithmic derivative of m. 
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To derive the critical exponent v we consider the loga- 
rithmic derivative of magnetization (see Fig. [6|) and use 

d In m i em — e-rn 



<9/3 



m 



(X ft em 



L 1 ^ + b Km - L yt+U 



(17) 

The equality in Eq. (|17[) implies that the equilibrium 
ensemble is well reproduced by our algorithm, which 
we expect to hold. Indeed, the obtained values for the 
exponent v are in excellent agreement with known ex- 
act, or best approximate results. The four-parameter fit 
to the form of the r.h.s. in (jTTJ) gave y T = 1.586(5), 
uj 4 = -2.0(2) and y T = 1.201(8), uj 4 = -1.0(3) for 3D 
Ising and 2D, q=3 case respectively. 
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Figure 6: (Color online) Logarithmic derivative of magneti- 
zation 



Table II: Corrections to scaling obtained by EIC algorithm 
(for reference see in the text) 





q = 2, 3D 
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-2.2(1) 
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-1.9(1) 
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-0.75(3) 


-0.91(4) 
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plotted versus the system size. 



I and II. It can be concluded that, with only a fraction of 
numerical effort needed for the SW algorithm, the critical 
behavior can be obtained very accurately and in excel- 
lent agreement with the best known or exact results. The 
location of the critical point can be obtained with rather 
high degree of accuracy as well. 

For better reliability of our results we included in our 
analysis the corrections to scaling. We observe that they 
differ substantially from the expected scaling corrections 
due to the irrelevant scaling fields and nonlinearities of 
scaling fields near criticality and should be attributed to 
the method itself. The leading irrelevant exponent for 2D 
q = 3 is exactly known to be yi = — |. For the 3D Ising 
case it is estimated to be j/j = —0.84(4) [22|. Our results 
show the scaling corrections that decay much faster than 
Hi in most cases, giving the exponents close to -2 or 
-1 in the 3D and 2D case respectively. The exception is 
the convergence of the Binder cumulant (W3 in Tab. II) 
with the correction exponent that can be compared to 
the value of leading irrelevant scaling field found by 
other authors. 



D. First-order phase transition 



Summary of critical exponents and corrections 
to scaling 



Table I: Leading critical parameters obtained by EIC (from 
L = 64 to 1024 for q = 3 2D and from L = 16 to 96 for 3D) 
compared to best known or exact results. 



q = 2, 3D 



= 3, 2D 



EIC best known" EIC 



exact* 



Pc{L 



oo)0.358097(l) 0.358098(3) 0.633975(1)0.6339745 . 
y T 1.586(5) 1.587(6) 1.201(8) f 
y h 2.4815(5) 2.4818(6) 1.8661(7) 1.86 
P/v 0.5180(4) 0.5182(6) 0.134(1) 0.13 
7/1/ 1.98(1) 1.964(1) 1.74(1) 1.73 



"Best known values from [22j and [2ll |. 

Critical exponents, and their leading corrections calcu- 
lated in preceding two sections are summarized in Tables 



The example for our study being the Potts model gives 
and opportunity to examine the applicability of our al- 
gorithm to the first order phase transitions. We make 
thus a small digression here, to show, without going into 
any details, that the same procedure, when applied to a 
first order phase transition gives a clear indication on its 
character. 

We first remind that for studies of the first-order tran- 
sitions by the IC algorithm the stopping rules based on 
the cluster mass are more appropriate Q and can be 
readily extended to EIC. Our intention here is to point 
out how the first-order phase transition manifests when 
the topological stopping rule is used. As the first-order 
transition is characterized by a finite hysteresis width, the 
construction of percolating clusters becomes possible for 
all the values of p{ a } within the hysteresis and thus the 
values of p^ sweep trough the entire width of the hys- 
teresis. This results with the plateau in the distribution 
w (P{<t}) who's width is not proportional to L~ 2 but tends 
to some constant value corresponding to the hysteresis 
width in the thermodynamic limit (see Fig. [7])- We illus- 
trate this on the two examples of the first-order transition 
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Figure 7: (Color online) Distributions io(p/ CT \) of P{ CT } pro- 
duced by EIC algorithm for the 3D Potts model in the case of 
strong (g = 5) and weak (q = 3) first-order phase transition. 



in the 3D Potts model. For q = 5, the first-order transi- 
tion is very strong and the figure clearly reveals the width 
of the hysteresis. For q = 3 which is characterized by a 
very weak first-order transition, larger lattices are needed 
to see the hysteresis. However, already for modest sizes, 
an inspection of the Figure [7| suggests that the width of 
the distribution will not vanish with some power-law in 
L, since its shape starts to qualitatively resemble the one 
for 3D, q = 5. 



IV. DYNAMICAL PROPERTIES OF EIC 
ALGORITHM 

The dynamics introduced by the EIC algorithm was 
studied by considering the autocorrelation function de- 
fined for an arbitrary observable O by 



r (*) 



O{t)O(0)-O 

&-o 2 



(18) 



The average is taken over the ensemble generated in the 
single run, after the thermalization, i.e. 



— - £ 0(t)0(t + t'), (19) 

(20) 



O(tf)O(0) 



t — t' 

max 1* 



t=l 



o» = - — ]T o(t) n , 



where t, and t' denote time in discrete integer units cor- 
responding to one MC step, and counting starts after 
equilibration. t max denotes the number of counted steps 
in the single run. The Eq. (|19p implies certain equi- 
librium properties of macroscopic variables such as the 
translational invariance in time and well defined temper- 
ature (or p). The self-regulating procedure of EIC con- 
tains deviations from these properties, but as these occur 



within bounds which we expect to preserve the equilib- 
rium ensemble we find justified to attempt a study of this 
process as the equilibrium one. We skip here the analysis 
of the nonequilibrium details of the process, which would 
require a different averaging procedure. Instead, we ap- 
ply the Eq. ()19p and just comment the differences found 
in comparison to the form of autocorrelation functions 
obtained in standard MC approaches. 

In standard MC approaches one may usually distin- 
guish three regimes of behavior of the autocorrelation 
functions [23|. For short times they behave as a sum of 
exponential decays. For longer times they adopt a single 
exponential decay of the form To{t) oc e~~ with a corre- 
lation time r, while for very large times the correlations 
vanish and the statistical errors take over. 

In Fig. [5] are displayed the autocorrelations in the EIC 
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Figure 8: (Color online) Magnetic autocorrelation function 
r m (t) for EIC, SW, and IC algorithms for 2D, q=3. On the 
inset of the figure are shown the integrated autocorrelation 
functions for the same cases. 

algorithm compared to those in the SW and IC algo- 
rithms on the example of the magnetic autocorrelation 
function in the case 2D q = 3. While the IC autocorre- 
lations fall extremely fast, the EIC correlations display 
similar decay as in the SW algorithm, with the difference 
that they acquire a negative sign at times of the order 
of N a . The slower decay of To (tj ) reflects the fact that 
IC dynamics has been slowed down by the requirement 
that the values of P{ a ) remain during N a MC steps in 
the narrow interval of 2 ■ v around the average of 
The fact that F(tj) acquires negative values can be re- 
lated to the abrupt switch of the mean value p i after 
every N a steps as a part of the self-regulating procedure. 
EIC algorithm has a tendency to " overshoot" the critical 
temperature and find itself on the other side of the crit- 
ical point in the following step. While O represents the 
average of the magnetization taken over the entire run, 
the average of magnetization over each single N a block 
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is different (larger of smaller than O) and varies accord- 
ing to p t . This explains why in the expression for the 
correlation function, (0(t) — O) and (0(t + t') — O) can 
in average be of the opposite sign, when t and t + t' be- 
long to the neighboring blocks of N a steps and produce 
anticorrclations. 

The semi- log plot displayed in Figure [H] shows that the 
autocorrelation functions indeed match well the assumed 
exponential form. This leads us to the question of the 
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Figure 9: (Color online) Energy autocorrelation functions 
r e (t) of q = 3, 2D case for several lattice sizes in the ex- 
ponential region. Dashed lines represent best fits. 

correlation time and the corresponding dynamic expo- 
nent z, which should satisfy the scaling relation 



T oc L z 



(21) 



Derivation of the correlation time is performed usually 
in two equivalent ways: by fitting the exponential decay 
(exponential correlation time), or by taking the integral 
r = J Q To(t)dt (integrated correlation time). 

The calculation of the exponential correlation time ap- 
plies to the present problem in a straightforward way, as 
already illustrated in Fig ([9]). The correlation times de- 
rived from the slopes of such semi- log plots and the corre- 
sponding dynamic exponents are presented in Table III, 
denoted by to j0 and zo, a respectively. The calculation of 
an integrated correlation time is obviously not applicable 
in the present case since Tq (t) is not positive in the entire 
domain of integration. In the inset of Fig. [8] we sketch 
the partial sums Jo 

(*) = Et'=o r o(*') for the three al- 
gorithms. While in the SW case the sum Ioit) rapidly 
saturates, for the EIC algorithm these sums tend to van- 
ish at infinity and exhibit maxima at the point where the 
To (i) turns negative. We thus consider instead the max- 
ima of these partial sums. Under certain approximation, 
if the correlation time is much smaller than the time in- 
terval N a /2, they may be understood as characteristic 
times To,b- Namely, if r <C N a /2, the exponential decay 



of autocorrelation functions dies out before the regime of 
anticorrelations sets in. The maxima of partial sums will 
then approximate well the (integrated) correlation time, 
while successfully leaving out of consideration the regime 
of anticorrelations. Our analysis shows that these max- 
ima (see Fig. I10[) exhibit indeed a power-law scaling with 
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Figure 10: (Color online) Integrals of magnetization auto- 
correlation functions I m (tk) = T^t k C 
several lattice sizes. 



r m (tj) of 3D Ising for 



an exponent very close to zo, a - However, if the condition 
t <C N a /2 is not fulfilled, such an approximation is not 
justified and we shall discuss the consequences later on. 
Both To, a (L) and To,b(L) were calculated for the autocor- 
relations of energy and magnetization and are presented 
in log-log plots versus the system size L in Fig, [TT1 while 
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Figure 11: (Color online) Characteristic times of energy and 
magnetization autocorrelation functions for the systems con- 
sidered, determined by two approaches described in the text. 
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the corresponding scaling exponents can be compared in 
Table III. 



Table III: Dynamical exponents of EIC algorithm (with N a = 
100, v = iiT*) 





q = 2, 3D 


q = 3, 2D 


EIC: Z m! a 


0.41(5) 


0.38(2) 


EIC: z m .b 


0.38(3) 


0.29(2) 


EIC: z e , a 


0.42(3) 


0.42(2) 


EIC: z e , b 


0.47(3) 


0.36(2) 


SW z b 


0.46(3) 


0.49(1) 


a/v 


0.172(1)° 


2/5 



"Best known values 123] . 
''Results cited in [25y5| 



Dynamical exponents derived from the energy and 
from the magnetization coincide up to the error mar- 
gins. We also observe a good matching between expo- 
nents zo, a and zo,b in the 3D case, while some discrep- 
ancy exists between these exponents in 2D. The reason 
for that can be related to the choice of the auxiliary pa- 
rameters. While the correlation times are several times 
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Figure 12: (Color online) Behavior of the autocorrelation 
functions of magnetization for U around 7V a = 100 (indicated 
by the vertical line) in the two cases considered and for various 
lattice sizes. 

larger in the 2D case, all the calculations presented in 
Table III were performed with the same N a = 100. As 
it may be observed in Fig. [12j for the 3D Ising case this 
choice is large enough to fulfill the condition r <C N a /2, 
and allow the autocorrelations to become negligible be- 
fore the effect of anticorrelations appears. In 2D case 
this condition is clearly not fulfilled, especially for larger 
lattice sizes, and the values T m ^ and the corresponding 
exponent cannot be related to the correlation time. It is 



interesting to discuss these results in context of the argu- 
ment by Li and Sokal [24[ , which sets the lower bound on 
the dynamical exponent z in cluster algorithms. Relating 
the minimal autocorrelation times to the heat capacity 
of the system, they argue that the lower bound for the 
dynamical exponent is given by 

z > a/v. (22) 

The results for z\eic summarized in Tab. Ill are com- 
pared to the most recent values of z\$w [IHIjIHI, and to 
the exponent ratio a/v. In both considered cases the dy- 
namical exponent appears to be comparable, but lower 
than the value for the SW algorithm. We also observe 
that the values tq a give results for the dynamical expo- 
nent conform to the Li - Sokal bound up to the statisti- 
cal error. In the 2D case, where discrepancy between the 
exponents zo, a and zo,b was observed, the dynamic ex- 
ponent is just about crossing the Li-Sokal limit, which 
can be understood as a sign that the system is close 
to leaving the equilibrium and taking the crossover to 
IC behavior. It has to be mentioned however, that as 
long as the width of temperature fluctuations is properly 
constrained this only seems to increase the noise in our 
results. No systematic errors were observed in the 2D re- 
sults for N a = 100, while the resulting ensemble retains 
the correct scaling, proper to the canonical ensemble as 
illustrated in Figure Q] 

In view of the above discussion we can drive some gen- 
eral conclusions about the choice of the auxiliary param- 
eters N a and v, and their impact on the results. 



A. Parameter N a 

It was briefly mentioned in the end of Sect. [TT| that 
the number N a must be large enough to prevent the IC 
dynamics from taking over the behavior of P{ a } fluctua- 
tions. This is the only restriction to the parameter N a . 
As discussed earlier in this Section, this starts to happen 
when N a is comparable to the value of the autocorrela- 
tion time of the system in consideration. Thus N a should 
be chosen to be at least more than two times larger than 
the autocorrelation time for any lattice size considered 
for a given system. When this is not fulfilled, there is a 
clear indication of the crossover to the IC dynamics in the 
behavior of the autocorrelation functions, as observed al- 
ready in the 2D case with N a = 100. When N a is reduced 
even further, the autocorrelation function changes the 
sign more than once, meaning that more than two values 
of Pj are correlated, which clearly indicates the domina- 
tion of IC dynamics over the EIC dynamics (Fig. [TB")) . 
Indeed, in these cases one finds that \fvar(p^ a }) cx L~ b 
where b < | and the initial requirement for fluctuations 
of P{ a } is not fulfilled. 
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Figure 13: (Color online) Autocorrelation functions of mag 
netization in the case 2D, q — 3 for different values of N a . 



B. Parameter v 

The effects of varying the parameter v can also be ob- 
served in the behavior of autocorrelation functions. We 
shall illustrate this on two examples. In the first example 
v was augmented by keeping the the power-law decay of 
the form L _d / 2 , but increasing the constant of propor- 
tionality v beyond the value compatible with the width 
following from the binomial distribution of Eq. In 
the second one, the parameter v lies inside the regime 
imposed by the EIC constraint, but decreases with size 
faster than L~ d / 2 . 

The examples are taken in the 2D, q = 3 case for the 
lattice sizes from L = 64 to 640, and using statistics of 
10 6 MC steps, while N a was increased to N a = 250 in 
order to avoid its contribution to nonequilibrium effects. 
For the case 2D, q = 3, the width of the distribution 
of P{ a ] corresponding to the binomial distribution in Eq. 

* « 0.19-i- 1 . 



would be : 



Vpc(i pc) n J ~ 019 .£-i We compare 
the results for the three different values of v. (a) v = 
(b) v = jq-L- 1 and (c) v = 1-L~ 15 . The results 
for the dynamical exponents are shown in the Tab. IV. 



Table IV: Dynamical exponents for the case 2D, q — 3 de- 
pending on v (with N a = 250) 







v=±-L 1 


v = l-L~ 15 


EIC: Zm,a 




0.39(3) 


0.56(2) 


EIC: z m ,b 


0.20(3) 


0.34(4) 


0.53(2) 


EIC: Ze,a 




0.41(3) 


0.52(2) 


EIC: z e , b 


0.33(4) 


0.38(4) 


0.54(2) 



Let us first comment the reference case (b) correspond- 
ing to the value of v used throughout this paper. When 
comparing the results with those in Table IV, obtained 



with N a — 100, one may see that by increasing N a the ex- 
ponent zo, a approaches closer to the Li-Sokal limit, while 
the discrepancy between exponents zo, a and zoj, dimin- 
ishes. This confirms the earlier discussion in this Section 
and the expectation that an increasing of N a would in- 
deed correct the deviation from the equilibrium behavior. 

In the case (a) the width v has the same power-law 
decay, but the constant factor v was taken larger than 
the value used in this paper and also larger than the 
one of the binomial distribution in Eq. ([5]). In spite of 
the correct power-law dependence, the crossover to the 
IC regime for lattice sizes considered here is observed. 
In autocorrelation functions, the regime of single expo- 
nential decay shrinks drastically and consequently times 
To.a could not be defined. We may still consider the 
quantity ro.b, since it still exhibits a power-law behav- 
ior, although it may not be interpreted as a correlation 
time. The corresponding exponent zo.b is much smaller 
than the Li-Sokal limit for the dynamical exponent of the 
equilibrium cluster algorithms. 

In contrast to this, in the example (c) with v = 1-L~ 15 
the exponential regime of the autocorrelation functions 
is longer then for v = ■ L^ 1 and consequently allows 
a more accurate calculation of both To, a and To,b, giving 
the same dynamical exponent up to the error margins. 
It is interesting to notice that by reducing the width pa- 
rameter v, the dynamic exponent z has increased. This 
provides a tool to vary the dynamical exponent contin- 
uously, by varying the exponent of power-law decay in 
v. 



V. DISCUSSION AND CONCLUSION 

In this paper we presented more detailed study of re- 
cently proposed EIC algorithm, designed to simultane- 
ously drive the system to criticality without its prior 
knowledge and to provide the correct critical exponents. 
Extensive numerical calculations performed on two typ- 
ical examples of the Potts model in two and three di- 
mensions confirmed that the method efficiently produces 
results for both critical exponents and critical temper- 
ature, with considerable accuracy that attains up to 4 
digits for critical exponents and six digits in the case of 
critical temperature. We conclude that the two auxiliary 
parameters of the algorithm, N a and v, may be varied 
in the wide range without producing systematic errors, 
and that their influence on the results when exceeding 
this range can be controlled by examining the behavior 
of the autocorrelation functions. We have also estimated 
the dynamic exponent of the EIC algorithm, and shown 
that it can be tuned by varying the auxiliary parameters. 
It is found to be larger than the one of the standard IC 
algorithm and generally smaller than the exponent of the 
SW algorithm. We also notice that, when the dynami- 
cal exponent is lowered and reaches beyond the Li-Sokal 
limit, the crossover to the nonequilibrium IC algorithm 
can be observed. 
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We observe that both in the SW method and in the 
present EIC calculations, the dynamical exponent does 
not change significantly with dimensionality and stays 
in the range between 0.4 and 0.5 in both cases, while 
the the Li-Sokal limit differs by more than twice between 
the two cases. The fact that within the EIC algorithm 
the dynamical exponent z may be tuned by changing the 
auxiliary parameters provides the means to make further 
analysis in order to better explain the reasons, still not 
well understood (26[, why the equilibrium cluster algo- 
rithms like SW in some cases (e.g. for 3D Ising), have 
much larger dynamical exponent z than required by the 
Li-Sokal limit. 

Due to its capability of self-tuning to the critical point 
there is a wide range of possible applications of the EIC 
algorithm to phase transitions, where the critical tem- 
perature is not known in advance. One of the most ap- 
pealing is certainly for simulations of phase transitions 



in presence of quenched disorder, where it provides the 
means to calculate critical properties at finite size critical 
temperature for each individual configuration at feasible 
time, which helps to overcome the problems related to 
the lack of self-averaging [27j • 

The algorithm could also be generalized to other, or 
more complex problems such as the first order phase tran- 
sition only briefly mentioned here, or to tricritical points, 
for which the extension of the standard IC algorithm was 
already applied 
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